function qr_iteration()
% QR_ITERATION QR迭代法求特征值
% 
% 演示QR迭代法及其变形用于计算所有特征值

fprintf('=== QR迭代法 ===\n\n');

%% 1. 基本QR迭代
fprintf('1. 基本QR迭代\n');
fprintf('迭代公式: A_k = Q_k * R_k, A_{k+1} = R_k * Q_k\n\n');

% 测试矩阵
A = [4, 1, 0; 1, 3, 1; 0, 1, 2];
fprintf('矩阵 A =\n');
disp(A);

% 精确特征值
eigenvals_exact = sort(eig(A), 'descend');
fprintf('精确特征值: [%.6f, %.6f, %.6f]\n', eigenvals_exact);

% QR迭代
max_iter = 50;
A_k = A;
fprintf('\nQR迭代过程:\n');
fprintf('迭代    对角元素                    收敛指标\n');
fprintf('----    --------                    --------\n');

for k = 1:max_iter
    [Q_k, R_k] = qr(A_k);
    A_k = R_k * Q_k;
    
    % 计算收敛指标（下三角元素的范数）
    convergence_measure = norm(tril(A_k, -1), 'fro');
    
    if k <= 10 || mod(k, 10) == 0 || convergence_measure < 1e-10
        diag_elements = diag(A_k);
        fprintf('%4d    [%6.4f, %6.4f, %6.4f]    %.2e\n', ...
                k, diag_elements, convergence_measure);
    end
    
    if convergence_measure < 1e-12
        fprintf('收敛于迭代 %d\n', k);
        break;
    end
end

fprintf('最终近似特征值: [%.6f, %.6f, %.6f]\n', sort(diag(A_k), 'descend'));

%% 2. 带位移的QR迭代
fprintf('\n2. 带位移的QR迭代\n');
fprintf('加速收敛的位移策略\n\n');

% 使用相同的测试矩阵
A_shift = A;
fprintf('带位移QR迭代:\n');
fprintf('迭代    位移        对角元素                收敛指标\n');
fprintf('----    ----        --------                --------\n');

for k = 1:30
    % 选择位移（使用右下角元素）
    n = size(A_shift, 1);
    shift = A_shift(n, n);
    
    % 带位移的QR分解
    [Q_shift, R_shift] = qr(A_shift - shift * eye(n));
    A_shift = R_shift * Q_shift + shift * eye(n);
    
    convergence_shift = norm(tril(A_shift, -1), 'fro');
    
    if k <= 8 || mod(k, 5) == 0 || convergence_shift < 1e-10
        diag_shift = diag(A_shift);
        fprintf('%4d    %6.4f    [%6.4f, %6.4f, %6.4f]    %.2e\n', ...
                k, shift, diag_shift, convergence_shift);
    end
    
    if convergence_shift < 1e-12
        fprintf('带位移QR收敛于迭代 %d\n', k);
        break;
    end
end

%% 3. Wilkinson位移策略
fprintf('\n3. Wilkinson位移策略\n');
fprintf('更智能的位移选择\n\n');

A_wilkinson = A;
fprintf('Wilkinson位移QR迭代:\n');

for k = 1:20
    n = size(A_wilkinson, 1);
    
    % Wilkinson位移：2×2子矩阵的特征值中更接近A(n,n)的那个
    if n >= 2
        submatrix = A_wilkinson(n-1:n, n-1:n);
        sub_eigenvals = eig(submatrix);
        [~, idx] = min(abs(sub_eigenvals - A_wilkinson(n, n)));
        wilkinson_shift = sub_eigenvals(idx);
    else
        wilkinson_shift = A_wilkinson(n, n);
    end
    
    [Q_w, R_w] = qr(A_wilkinson - wilkinson_shift * eye(n));
    A_wilkinson = R_w * Q_w + wilkinson_shift * eye(n);
    
    convergence_w = norm(tril(A_wilkinson, -1), 'fro');
    
    if k <= 6 || convergence_w < 1e-10
        fprintf('迭代 %2d: 位移 = %8.6f, 收敛指标 = %.2e\n', ...
                k, wilkinson_shift, convergence_w);
    end
    
    if convergence_w < 1e-12
        fprintf('Wilkinson位移QR收敛于迭代 %d\n', k);
        break;
    end
end

%% 4. 双重位移QR迭代
fprintf('\n4. 双重位移QR迭代\n');
fprintf('处理复特征值的实算法\n\n');

% 构造有复特征值的矩阵
A_complex = [0, 1, 0; -1, 0, 0; 0, 0, 2];
fprintf('有复特征值的矩阵:\n');
disp(A_complex);

eigenvals_complex = eig(A_complex);
fprintf('精确特征值: %.4f±%.4fi, %.4f\n', ...
        real(eigenvals_complex(1)), imag(eigenvals_complex(1)), eigenvals_complex(3));

% 双重位移QR（简化版本）
A_double = A_complex;
fprintf('双重位移QR迭代:\n');

for k = 1:15
    n = size(A_double, 1);
    
    if n >= 3
        % 使用右下角2×2子矩阵的特征值作为双重位移
        submat = A_double(n-1:n, n-1:n);
        shifts = eig(submat);
        
        % 应用双重位移（简化实现）
        if imag(shifts(1)) ~= 0
            % 复共轭位移对
            sigma1 = shifts(1);
            sigma2 = conj(sigma1);
            
            % (A - σ₁I)(A - σ₂I) = A² - (σ₁+σ₂)A + σ₁σ₂I
            s = sigma1 + sigma2;  % 实数
            t = sigma1 * sigma2;  % 实数
            
            M = A_double^2 - s*A_double + t*eye(n);
            [Q_double, R_double] = qr(M);
            A_double = Q_double' * A_double * Q_double;
        else
            % 实位移
            [Q_double, R_double] = qr(A_double - shifts(1)*eye(n));
            A_double = R_double * Q_double + shifts(1)*eye(n);
        end
    else
        [Q_double, R_double] = qr(A_double);
        A_double = R_double * Q_double;
    end
    
    conv_double = norm(tril(A_double, -1), 'fro');
    
    if k <= 8 || conv_double < 1e-8
        fprintf('迭代 %2d: 收敛指标 = %.2e\n', k, conv_double);
    end
    
    if conv_double < 1e-10
        break;
    end
end

%% 5. Hessenberg预处理
fprintf('\n5. Hessenberg预处理\n');
fprintf('将矩阵化为上Hessenberg形式加速QR迭代\n\n');

% 较大的测试矩阵
A_large = randn(6, 6);
A_large = (A_large + A_large') / 2;  % 对称化
fprintf('测试矩阵 (%d×%d)\n', size(A_large));

% 直接QR迭代
tic;
A_direct = A_large;
for k = 1:100
    [Q_d, R_d] = qr(A_direct);
    A_direct = R_d * Q_d;
    if norm(tril(A_direct, -1), 'fro') < 1e-10
        break;
    end
end
time_direct = toc;
fprintf('直接QR迭代: %d次迭代, %.4f秒\n', k, time_direct);

% Hessenberg + QR迭代
tic;
[P, H] = hess(A_large);  % 化为Hessenberg形式
A_hess = H;
for k = 1:100
    [Q_h, R_h] = qr(A_hess);
    A_hess = R_h * Q_h;
    if norm(tril(A_hess, -2), 'fro') < 1e-10  % Hessenberg矩阵只需检查下下三角
        break;
    end
end
time_hess = toc;
fprintf('Hessenberg+QR: %d次迭代, %.4f秒\n', k, time_hess);
fprintf('加速比: %.2f\n', time_direct / time_hess);

%% 6. 收缩技术
fprintf('\n6. 收缩技术\n');
fprintf('逐个分离特征值\n\n');

A_deflate = [5, 1, 0, 0; 1, 4, 1, 0; 0, 1, 3, 1; 0, 0, 1, 2];
fprintf('测试矩阵 A (%d×%d):\n', size(A_deflate));
disp(A_deflate);

eigenvals_deflate = sort(eig(A_deflate), 'descend');
fprintf('精确特征值: [%.4f, %.4f, %.4f, %.4f]\n', eigenvals_deflate);

% 模拟收缩过程
A_current = A_deflate;
found_eigenvals = [];

fprintf('收缩过程:\n');
for deflation_step = 1:4
    n_current = size(A_current, 1);
    fprintf('步骤 %d: 处理 %d×%d 矩阵\n', deflation_step, n_current, n_current);
    
    if n_current == 1
        found_eigenvals = [found_eigenvals; A_current(1, 1)];
        fprintf('找到特征值: %.6f\n', A_current(1, 1));
        break;
    end
    
    % QR迭代直到收敛
    for iter = 1:50
        [Q_def, R_def] = qr(A_current);
        A_current = R_def * Q_def;
        
        % 检查是否可以收缩
        if abs(A_current(n_current, n_current-1)) < 1e-10
            % 分离出一个特征值
            found_eigenvals = [found_eigenvals; A_current(n_current, n_current)];
            fprintf('找到特征值: %.6f (迭代 %d 次)\n', ...
                    A_current(n_current, n_current), iter);
            
            % 收缩矩阵
            A_current = A_current(1:n_current-1, 1:n_current-1);
            break;
        end
    end
end

fprintf('找到的特征值: [%.4f, %.4f, %.4f, %.4f]\n', sort(found_eigenvals, 'descend'));

%% 7. 实对称矩阵的QR迭代
fprintf('\n7. 实对称矩阵的QR迭代\n');
fprintf('利用对称性的特殊算法\n\n');

% 对称矩阵
A_sym = [4, 1, 0; 1, 3, 1; 0, 1, 2];
fprintf('对称矩阵 A =\n');
disp(A_sym);

% 三对角化
[V_tri, T_tri] = schur(A_sym);  % 使用Schur分解得到三对角形式
fprintf('三对角化后的矩阵:\n');
disp(T_tri);

% 对称三对角矩阵的QR迭代收敛更快
A_tri = T_tri;
fprintf('对称三对角QR迭代:\n');

for k = 1:20
    [Q_tri, R_tri] = qr(A_tri);
    A_tri = R_tri * Q_tri;
    
    off_diag_norm = norm(A_tri - diag(diag(A_tri)), 'fro');
    
    if k <= 8 || off_diag_norm < 1e-10
        fprintf('迭代 %2d: 非对角元素范数 = %.2e\n', k, off_diag_norm);
    end
    
    if off_diag_norm < 1e-12
        fprintf('对称QR收敛于迭代 %d\n', k);
        break;
    end
end

%% 8. 性能比较
fprintf('\n8. 性能比较\n');
fprintf('不同QR迭代策略的性能对比\n\n');

test_sizes = [5, 10, 20];
fprintf('矩阵规模    基本QR    带位移QR    Hessenberg+QR\n');
fprintf('--------    ------    --------    -------------\n');

for i = 1:length(test_sizes)
    n_test = test_sizes(i);
    A_test = randn(n_test, n_test);
    
    % 基本QR迭代
    tic;
    A_basic = A_test;
    for k = 1:200
        [Q_basic, R_basic] = qr(A_basic);
        A_basic = R_basic * Q_basic;
        if norm(tril(A_basic, -1), 'fro') < 1e-8
            break;
        end
    end
    time_basic = toc;
    iter_basic = k;
    
    % 带位移QR迭代
    tic;
    A_shifted = A_test;
    for k = 1:200
        shift_test = A_shifted(end, end);
        [Q_shifted, R_shifted] = qr(A_shifted - shift_test*eye(n_test));
        A_shifted = R_shifted * Q_shifted + shift_test*eye(n_test);
        if norm(tril(A_shifted, -1), 'fro') < 1e-8
            break;
        end
    end
    time_shifted = toc;
    iter_shifted = k;
    
    % Hessenberg + QR
    tic;
    [P_test, H_test] = hess(A_test);
    A_hess_test = H_test;
    for k = 1:200
        [Q_hess_test, R_hess_test] = qr(A_hess_test);
        A_hess_test = R_hess_test * Q_hess_test;
        if norm(tril(A_hess_test, -2), 'fro') < 1e-8
            break;
        end
    end
    time_hess_test = toc;
    iter_hess = k;
    
    fprintf('%6d      %6.4f    %8.4f    %13.4f\n', ...
            n_test, time_basic, time_shifted, time_hess_test);
    fprintf('        (%3d次)    (%3d次)      (%3d次)\n', ...
            iter_basic, iter_shifted, iter_hess);
end

%% 9. 应用：振动分析
fprintf('\n9. 应用：振动分析\n');
fprintf('使用QR迭代求解振动系统的固有频率\n\n');

% 质量和刚度矩阵
M = diag([2, 1, 2]);
K = [3, -1, 0; -1, 2, -1; 0, -1, 1];

fprintf('质量矩阵 M =\n'); disp(M);
fprintf('刚度矩阵 K =\n'); disp(K);

% 广义特征值问题转换为标准特征值问题
% K*v = λ*M*v => M^(-1)*K*v = λ*v
A_vibration = M \ K;
fprintf('标准特征值问题矩阵 M^(-1)*K =\n');
disp(A_vibration);

% 使用QR迭代求解
A_vib = A_vibration;
fprintf('振动系统QR迭代:\n');

for k = 1:30
    [Q_vib, R_vib] = qr(A_vib);
    A_vib = R_vib * Q_vib;
    
    if norm(tril(A_vib, -1), 'fro') < 1e-10
        break;
    end
end

frequencies = sqrt(sort(diag(A_vib), 'descend'));
fprintf('固有频率 (rad/s): [%.4f, %.4f, %.4f]\n', frequencies);

% 与精确解比较
exact_frequencies = sqrt(sort(eig(A_vibration), 'descend'));
fprintf('精确固有频率: [%.4f, %.4f, %.4f]\n', exact_frequencies);
fprintf('误差: [%.2e, %.2e, %.2e]\n', abs(frequencies - exact_frequencies));

end